# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import numpy as np
import functools
from hysop.tools.htypes import check_instance, first_not_None
from hysop.tools.decorators import debug
from hysop.tools.numpywrappers import npw
from hysop.backend.device.opencl.opencl_operator import OpenClOperator
from hysop.core.graph.graph import op_apply
from hysop.fields.continuous_field import Field
from hysop.parameters.parameter import Parameter
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
from hysop.operator.base.spatial_filtering import (
RemeshRestrictionFilterBase,
SpectralRestrictionFilterBase,
SubgridRestrictionFilterBase,
PolynomialInterpolationFilterBase,
PolynomialRestrictionFilterBase,
)
from hysop.backend.device.opencl.opencl_symbolic import OpenClSymbolic
from hysop.backend.device.opencl.opencl_copy_kernel_launchers import (
OpenClCopyBufferRectLauncher,
)
from hysop.backend.device.opencl.opencl_kernel_launcher import OpenClKernelListLauncher
from hysop.backend.device.opencl.opencl_elementwise import (
OpenClElementwiseKernelGenerator,
)
from hysop.symbolic import local_indices_symbols
from hysop.symbolic.relational import Assignment
[docs]
class OpenClPolynomialInterpolationFilter(
PolynomialInterpolationFilterBase, OpenClOperator
):
[docs]
@debug
def discretize(self):
if self.discretized:
return
super().discretize()
dFin = self.dFin
dFout = self.dFout
gr = self.grid_ratio
dim = dFin.dim
assert dFin.is_scalar
assert dFout.is_scalar
assert self.subgrid_interpolator.gr == gr
ekg = self.elementwise_kernel_generator
Wr = self.subgrid_interpolator.Wr
n = self.subgrid_interpolator.n
ghosts = np.asarray(self.subgrid_interpolator.ghosts)
I = np.asarray(local_indices_symbols[:dim][::-1])
fin, fout = ekg.dfields_to_ndbuffers(dFin, dFout)
Fin = ekg.symbolic_tmp_scalars("F", shape=n, dtype=dFin.dtype)
Fout_values = Wr.dot(Fin.ravel()).reshape(gr)
exprs = ()
for idx in np.ndindex(*n):
e = Assignment(Fin[idx], fin(I + idx - ghosts))
exprs += (e,)
for idx in np.ndindex(*gr):
e = Assignment(fout(gr * I + idx), Fout_values[idx])
exprs += (e,)
kname = f"interpolate_grid_{self.polynomial_interpolation_method}".lower()
interpolate_grid_kernel, _ = ekg.elementwise_kernel(
kname, *exprs, compute_resolution=self.iter_shape, debug=False
)
exchange_ghosts = self.dFout.exchange_ghosts(build_launcher=True)
kl = OpenClKernelListLauncher(name=kname, profiler=self._profiler)
kl += interpolate_grid_kernel
kl += exchange_ghosts
self.execute_kernels = functools.partial(kl, queue=self.cl_env.default_queue)
@op_apply
def apply(self, **kwds):
super().apply(**kwds)
evt = self.execute_kernels()
[docs]
class OpenClPolynomialRestrictionFilter(
PolynomialRestrictionFilterBase, OpenClOperator
):
[docs]
@debug
def discretize(self):
if self.discretized:
return
super().discretize()
dFin = self.dFin
dFout = self.dFout
gr = self.grid_ratio
dim = dFin.dim
assert dFin.is_scalar
assert dFout.is_scalar
assert self.subgrid_restrictor.gr == gr
ekg = self.elementwise_kernel_generator
Rr = self.subgrid_restrictor.Rr / self.subgrid_restrictor.gr
ghosts = np.asarray(self.subgrid_restrictor.ghosts)
I = np.asarray(local_indices_symbols[:dim][::-1])
fin, fout = ekg.dfields_to_ndbuffers(dFin, dFout)
def gen_inputs(*idx):
return fin(gr * I + idx - ghosts)
input_values = np.asarray(
tuple(map(gen_inputs, np.ndindex(*Rr.shape)))
).reshape(Rr.shape)
output_value = (Rr * input_values).sum()
e = Assignment(fout(I), output_value)
exprs = (e,)
kname = f"restrict_grid_{self.polynomial_interpolation_method}".lower()
restriction_grid_kernel, _ = ekg.elementwise_kernel(
kname, *exprs, compute_resolution=self.iter_shape, debug=False
)
exchange_ghosts = self.dFout.exchange_ghosts(build_launcher=True)
kl = OpenClKernelListLauncher(name=kname, profiler=self._profiler)
kl += restriction_grid_kernel
kl += exchange_ghosts
self.execute_kernels = functools.partial(kl, queue=self.cl_env.default_queue)
@op_apply
def apply(self, **kwds):
super().apply(**kwds)
evt = self.execute_kernels()
[docs]
class OpenClSubgridRestrictionFilter(SubgridRestrictionFilterBase, OpenClSymbolic):
"""
OpenCL implementation for lowpass spatial filtering: small grid -> coarse grid
using the subgrid method.
"""
def __init__(self, **kwds):
super().__init__(**kwds)
Fin = self.Fin
Fout = self.Fout
dim = Fin.dim
assert Fin.is_scalar
assert Fout.is_scalar
# We do not know the grid ratio and array strides before discretization.
# so we defer the initialization of those integers with symbolic constants.
(symbolic_input_buffer,) = self.symbolic_buffers("fine_grid")
symbolic_output_buffer = self.Fout.s()
symbolic_grid_ratio = self.symbolic_constants("gr", count=dim, dtype=npw.int32)
symbolic_input_strides = self.symbolic_constants(
"is", count=dim, dtype=npw.int32
)
symbolic_input_ghosts = self.symbolic_constants(
"gs", count=dim, dtype=npw.int32
)
I = local_indices_symbols[:dim][::-1]
read_idx = npw.dot(
symbolic_input_strides,
npw.add(npw.multiply(symbolic_grid_ratio, I), symbolic_input_ghosts),
)
expr = Assignment(symbolic_output_buffer, symbolic_input_buffer[read_idx])
self.require_symbolic_kernel("extract_subgrid", expr)
self.symbolic_input_buffer = symbolic_input_buffer
self.symbolic_output_buffer = symbolic_output_buffer
self.symbolic_grid_ratio = symbolic_grid_ratio
self.symbolic_input_strides = symbolic_input_strides
self.symbolic_input_ghosts = symbolic_input_ghosts
[docs]
@debug
def setup(self, work):
dFin, dFout = self.dFin, self.dFout
ibuffer, obuffer = dFin.sbuffer, dFout.sbuffer
self.symbolic_input_buffer.bind_memory_object(ibuffer)
for i in range(dFin.dim):
self.symbolic_grid_ratio[i].bind_value(self.grid_ratio[i])
self.symbolic_input_strides[i].bind_value(
ibuffer.strides[i] // ibuffer.dtype.itemsize
)
self.symbolic_input_ghosts[i].bind_value(dFin.ghosts[i])
super().setup(work)
(extract_subgrid, _) = self.symbolic_kernels["extract_subgrid"]
exchange_ghosts = self.dFout.exchange_ghosts(build_launcher=True)
kl = OpenClKernelListLauncher(name="extract_subgrid")
kl += extract_subgrid
kl += exchange_ghosts
self.execute_kernels = functools.partial(kl, queue=self.cl_env.default_queue)
@op_apply
def apply(self, **kwds):
evt = self.execute_kernels()
[docs]
class OpenClSpectralRestrictionFilter(SpectralRestrictionFilterBase, OpenClOperator):
"""
OpenCL implementation for lowpass spatial filtering: small grid -> coarse grid
using the spectral method.
"""
def _compute_scaling_coefficient(self):
kernel_generator = OpenClElementwiseKernelGenerator(
cl_env=self.cl_env, kernel_config=self.kernel_config
)
# Kernels to copy src_slices to dst_slices (windowing operation on frequencies)
kl = OpenClKernelListLauncher(name="lowpass_filter")
for src_slc, dst_slc in zip(*self.fslices):
kl += OpenClCopyBufferRectLauncher.from_slices(
"copy", src=self.FIN[src_slc], dst=self.FOUT[dst_slc]
)
# Now we compute the scaling coefficient of the filter
# self.Ft.input_buffer is just a pypencl.Array so we need to use
# the kernel_generator to fill ones and use explicit copy kernels
# This seems to be the only solution to fill a non C-contiguous
# OpenClinput_buffer with ones.
buf = self.Ft.input_buffer
(buf,) = kernel_generator.arrays_to_symbols(buf)
expr = Assignment(buf, 1)
fill_ones, _ = kernel_generator.elementwise_kernel("fill_ones", expr)
fill_ones(queue=self.cl_env.default_queue)
self.Ft(simulation=False)
kl(queue=self.cl_env.default_queue) # here we apply unscaled filter
self.Bt(simulation=False)
# Here we get the coefficient
scaling = 1.0 / self.Bt.output_buffer[(0,) * self.FOUT.ndim].get()
# Now we can finally build the filter scaling kernel
(fout,) = kernel_generator.arrays_to_symbols(self.FOUT)
expr = Assignment(fout, scaling * fout)
scale, _ = kernel_generator.elementwise_kernel("scale", expr)
kl += scale
# we store the filtering kernel list for the setup step
self.filter = functools.partial(kl, queue=self.cl_env.default_queue)
# finally build ghost exchanger
exchange_ghosts = self.dFout.exchange_ghosts(build_launcher=True)
if exchange_ghosts is not None:
exchange_ghosts = functools.partial(
exchange_ghosts, queue=self.cl_env.default_queue
)
self.exchange_ghosts = exchange_ghosts
return scaling
@op_apply
def apply(self, **kwds):
"""Apply spectral filter (which is just a square window centered on low frequencies)."""
super().apply(**kwds)
evt = self.Ft(**kwds)
evt = self.filter()
evt = self.Bt(**kwds)
if self.exchange_ghosts is not None:
evt = self.exchange_ghosts()